Whats in V3: GAM models using Humidity, Temperature, PM2.5 and NPIs

Climate variables (temperature, humidity, PM2.5) lagged so that daily cases are predicted using climate variables from n=climate_lag days prior

NPI lagged so that daily cases are predicted using NPI data from n=0,…,num_iter days prior

the climate splines should be more or less the same since we’re keeping climate_lag constant i think what we’re looking for is an appropriate lag for npi which should change a lot over the num_iter GAMs

If you run this you’ll see that infections decrease with temperature and increase with PM2.5 Humidity is kind of a wildcard. This is overall true for all the cities (florence is a little weird) and for all the values of climate lag I tried (3,4,5,6)

We start by loading in our data tables

require(mgcv)
bologna <- read.csv("bologna_full.csv")
bologna <- na.omit(bologna)

modena <- read.csv("modena_full.csv")
modena <- na.omit(modena)

florence <- read.csv("florence_full.csv")
florence <- na.omit(florence)

rome <- read.csv("rome_full.csv")
rome <- na.omit(rome)

num_iter <- 30
climate_lag <- 6

Bologna:

for(i in 0:num_iter) {
num_lag <- i
bologna <- read.csv("bologna_full.csv") #Read in file
#Push temperature&humidity data back climate_lag days by adding empty data in front
#rep(NA, n) creates a vector of length n filled with NA
bologna$Temperature.1 <- c(rep(NA, climate_lag),  bologna$Temperature[(1):(nrow(bologna)-climate_lag)]) 
bologna$Humidity.1 <- c(rep(NA, climate_lag), bologna$Humidity[(1):(nrow(bologna)-climate_lag)])
bologna$pm25.1 <- c(rep(NA, climate_lag), bologna$pm25[(1):(nrow(bologna)-climate_lag)])


#push npi back num_lag days
bologna$npi.1 <- c(rep(NA, num_lag), bologna$npi[(1):(nrow(bologna)-num_lag)])
#Omit the new NAs
bologna <- na.omit(bologna)

infections.gam <- gam(Daily ~ s(Temperature.1) + s(Humidity.1) + s(pm25.1) + s(npi.1, k=4), family=poisson(link="log"), data=bologna)
plot(infections.gam,scale=0,se=2, shade=TRUE,pages=1)
title(main = i)
}

Modena:

for(i in 0:num_iter) {
num_lag <- i
modena <- read.csv("modena_full.csv") #Read in file
#Push temperature&humidity data back climate_lag days by adding empty data in front
#rep(NA, n) creates a vector of length n filled with NA
modena$Temperature.1 <- c(rep(NA, climate_lag),  modena$Temperature[(1):(nrow(modena)-climate_lag)]) 
modena$Humidity.1 <- c(rep(NA, climate_lag), modena$Humidity[(1):(nrow(modena)-climate_lag)])
modena$pm25.1 <- c(rep(NA, climate_lag), modena$pm25[(1):(nrow(modena)-climate_lag)])

#push npi back num_lag days
modena$npi.1 <- c(rep(NA, num_lag), modena$npi[(1):(nrow(modena)-num_lag)])
#Omit the new NAs
modena <- na.omit(modena)

infections.gam <- gam(Daily ~ s(Temperature.1) + s(Humidity.1) + s(pm25.1) + s(npi.1, k=4), family=poisson(link="log"), data=modena)
plot(infections.gam,scale=0,se=2, shade=TRUE,pages=1)
title(main = i)
}

NA

Florence:

for(i in 0:num_iter) {
num_lag <- i
florence <- read.csv("florence_full.csv") #Read in file
#Push temperature&humidity data back climate_lag days by adding empty data in front
#rep(NA, n) creates a vector of length n filled with NA
florence$Temperature.1 <- c(rep(NA, climate_lag),  florence$Temperature[(1):(nrow(florence)-climate_lag)]) 
florence$Humidity.1 <- c(rep(NA, climate_lag), florence$Humidity[(1):(nrow(florence)-climate_lag)])
florence$pm25.1 <- c(rep(NA, climate_lag), florence$pm25[(1):(nrow(florence)-climate_lag)])

#push npi back num_lag days
florence$npi.1 <- c(rep(NA, num_lag), florence$npi[(1):(nrow(florence)-num_lag)])
#Omit the new NAs
florence <- na.omit(florence)

infections.gam <- gam(Daily ~ s(Temperature.1) + s(Humidity.1) + s(pm25.1) + s(npi.1, k=4), family=poisson(link="log"), data=florence)
plot(infections.gam,scale=0,se=2, shade=TRUE,pages=1)
title(main = i)
}

Rome:

for(i in 0:num_iter) {
num_lag <- i
rome <- read.csv("rome_full.csv") #Read in file
#Push temperature&humidity data back climate_lag days by adding empty data in front
#rep(NA, n) creates a vector of length n filled with NA
rome$Temperature.1 <- c(rep(NA, climate_lag),  rome$Temperature[(1):(nrow(rome)-climate_lag)]) 
rome$Humidity.1 <- c(rep(NA, climate_lag), rome$Humidity[(1):(nrow(rome)-climate_lag)])
rome$pm25.1 <- c(rep(NA, climate_lag), rome$pm25[(1):(nrow(rome)-climate_lag)])

#push npi back num_lag days
rome$npi.1 <- c(rep(NA, num_lag), rome$npi[(1):(nrow(rome)-num_lag)])
#Omit the new NAs
rome <- na.omit(rome)

infections.gam <- gam(Daily ~ s(Temperature.1) + s(Humidity.1) + s(pm25.1)+s(npi.1, k=4), family=poisson(link="log"), data=rome)
plot(infections.gam,scale=0,se=2, shade=TRUE,pages=1)
title(main = i)
}

LS0tCnRpdGxlOiAiUiBOb3RlYm9vayIKb3V0cHV0OiBodG1sX25vdGVib29rCi0tLQpXaGF0cyBpbiBWMzoKR0FNIG1vZGVscyB1c2luZyBIdW1pZGl0eSwgVGVtcGVyYXR1cmUsIFBNMi41IGFuZCBOUElzCgpDbGltYXRlIHZhcmlhYmxlcyAodGVtcGVyYXR1cmUsIGh1bWlkaXR5LCBQTTIuNSkgbGFnZ2VkIHNvIHRoYXQgZGFpbHkgY2FzZXMgYXJlIHByZWRpY3RlZCB1c2luZyBjbGltYXRlIHZhcmlhYmxlcyBmcm9tIG49Y2xpbWF0ZV9sYWcgZGF5cyBwcmlvcgoKTlBJIGxhZ2dlZCBzbyB0aGF0IGRhaWx5IGNhc2VzIGFyZSBwcmVkaWN0ZWQgdXNpbmcgTlBJIGRhdGEgZnJvbSBuPTAsLi4uLG51bV9pdGVyIGRheXMgcHJpb3IKCnRoZSBjbGltYXRlIHNwbGluZXMgc2hvdWxkIGJlIG1vcmUgb3IgbGVzcyB0aGUgc2FtZSBzaW5jZSB3ZSdyZSBrZWVwaW5nIGNsaW1hdGVfbGFnIGNvbnN0YW50CmkgdGhpbmsgd2hhdCB3ZSdyZSBsb29raW5nIGZvciBpcyBhbiBhcHByb3ByaWF0ZSBsYWcgZm9yIG5waSB3aGljaCBzaG91bGQgY2hhbmdlIGEgbG90IG92ZXIgdGhlIG51bV9pdGVyIEdBTXMKCklmIHlvdSBydW4gdGhpcyB5b3UnbGwgc2VlIHRoYXQgaW5mZWN0aW9ucyBkZWNyZWFzZSB3aXRoIHRlbXBlcmF0dXJlIGFuZCBpbmNyZWFzZSB3aXRoIFBNMi41Ckh1bWlkaXR5IGlzIGtpbmQgb2YgYSB3aWxkY2FyZC4gClRoaXMgaXMgb3ZlcmFsbCB0cnVlIGZvciBhbGwgdGhlIGNpdGllcyAoZmxvcmVuY2UgaXMgYSBsaXR0bGUgd2VpcmQpIGFuZCBmb3IgYWxsIHRoZSB2YWx1ZXMgb2YgY2xpbWF0ZSBsYWcgSSB0cmllZCAoMyw0LDUsNikKCgpXZSBzdGFydCBieSBsb2FkaW5nIGluIG91ciBkYXRhIHRhYmxlcwoKYGBge3J9CnJlcXVpcmUobWdjdikKYm9sb2duYSA8LSByZWFkLmNzdigiYm9sb2duYV9mdWxsLmNzdiIpCmJvbG9nbmEgPC0gbmEub21pdChib2xvZ25hKQoKbW9kZW5hIDwtIHJlYWQuY3N2KCJtb2RlbmFfZnVsbC5jc3YiKQptb2RlbmEgPC0gbmEub21pdChtb2RlbmEpCgpmbG9yZW5jZSA8LSByZWFkLmNzdigiZmxvcmVuY2VfZnVsbC5jc3YiKQpmbG9yZW5jZSA8LSBuYS5vbWl0KGZsb3JlbmNlKQoKcm9tZSA8LSByZWFkLmNzdigicm9tZV9mdWxsLmNzdiIpCnJvbWUgPC0gbmEub21pdChyb21lKQoKbnVtX2l0ZXIgPC0gMzAKY2xpbWF0ZV9sYWcgPC0gNgpgYGAKCgpCb2xvZ25hOgpgYGB7cn0KZm9yKGkgaW4gMDpudW1faXRlcikgewpudW1fbGFnIDwtIGkKYm9sb2duYSA8LSByZWFkLmNzdigiYm9sb2duYV9mdWxsLmNzdiIpICNSZWFkIGluIGZpbGUKI1B1c2ggdGVtcGVyYXR1cmUmaHVtaWRpdHkgZGF0YSBiYWNrIGNsaW1hdGVfbGFnIGRheXMgYnkgYWRkaW5nIGVtcHR5IGRhdGEgaW4gZnJvbnQKI3JlcChOQSwgbikgY3JlYXRlcyBhIHZlY3RvciBvZiBsZW5ndGggbiBmaWxsZWQgd2l0aCBOQQpib2xvZ25hJFRlbXBlcmF0dXJlLjEgPC0gYyhyZXAoTkEsIGNsaW1hdGVfbGFnKSwgIGJvbG9nbmEkVGVtcGVyYXR1cmVbKDEpOihucm93KGJvbG9nbmEpLWNsaW1hdGVfbGFnKV0pIApib2xvZ25hJEh1bWlkaXR5LjEgPC0gYyhyZXAoTkEsIGNsaW1hdGVfbGFnKSwgYm9sb2duYSRIdW1pZGl0eVsoMSk6KG5yb3coYm9sb2duYSktY2xpbWF0ZV9sYWcpXSkKYm9sb2duYSRwbTI1LjEgPC0gYyhyZXAoTkEsIGNsaW1hdGVfbGFnKSwgYm9sb2duYSRwbTI1WygxKToobnJvdyhib2xvZ25hKS1jbGltYXRlX2xhZyldKQoKCiNwdXNoIG5waSBiYWNrIG51bV9sYWcgZGF5cwpib2xvZ25hJG5waS4xIDwtIGMocmVwKE5BLCBudW1fbGFnKSwgYm9sb2duYSRucGlbKDEpOihucm93KGJvbG9nbmEpLW51bV9sYWcpXSkKI09taXQgdGhlIG5ldyBOQXMKYm9sb2duYSA8LSBuYS5vbWl0KGJvbG9nbmEpCgppbmZlY3Rpb25zLmdhbSA8LSBnYW0oRGFpbHkgfiBzKFRlbXBlcmF0dXJlLjEpICsgcyhIdW1pZGl0eS4xKSArIHMocG0yNS4xKSArIHMobnBpLjEsIGs9NCksIGZhbWlseT1wb2lzc29uKGxpbms9ImxvZyIpLCBkYXRhPWJvbG9nbmEpCnBsb3QoaW5mZWN0aW9ucy5nYW0sc2NhbGU9MCxzZT0yLCBzaGFkZT1UUlVFLHBhZ2VzPTEpCnRpdGxlKG1haW4gPSBpKQp9CmBgYAoKTW9kZW5hOgpgYGB7cn0KZm9yKGkgaW4gMDpudW1faXRlcikgewpudW1fbGFnIDwtIGkKbW9kZW5hIDwtIHJlYWQuY3N2KCJtb2RlbmFfZnVsbC5jc3YiKSAjUmVhZCBpbiBmaWxlCiNQdXNoIHRlbXBlcmF0dXJlJmh1bWlkaXR5IGRhdGEgYmFjayBjbGltYXRlX2xhZyBkYXlzIGJ5IGFkZGluZyBlbXB0eSBkYXRhIGluIGZyb250CiNyZXAoTkEsIG4pIGNyZWF0ZXMgYSB2ZWN0b3Igb2YgbGVuZ3RoIG4gZmlsbGVkIHdpdGggTkEKbW9kZW5hJFRlbXBlcmF0dXJlLjEgPC0gYyhyZXAoTkEsIGNsaW1hdGVfbGFnKSwgIG1vZGVuYSRUZW1wZXJhdHVyZVsoMSk6KG5yb3cobW9kZW5hKS1jbGltYXRlX2xhZyldKSAKbW9kZW5hJEh1bWlkaXR5LjEgPC0gYyhyZXAoTkEsIGNsaW1hdGVfbGFnKSwgbW9kZW5hJEh1bWlkaXR5WygxKToobnJvdyhtb2RlbmEpLWNsaW1hdGVfbGFnKV0pCm1vZGVuYSRwbTI1LjEgPC0gYyhyZXAoTkEsIGNsaW1hdGVfbGFnKSwgbW9kZW5hJHBtMjVbKDEpOihucm93KG1vZGVuYSktY2xpbWF0ZV9sYWcpXSkKCiNwdXNoIG5waSBiYWNrIG51bV9sYWcgZGF5cwptb2RlbmEkbnBpLjEgPC0gYyhyZXAoTkEsIG51bV9sYWcpLCBtb2RlbmEkbnBpWygxKToobnJvdyhtb2RlbmEpLW51bV9sYWcpXSkKI09taXQgdGhlIG5ldyBOQXMKbW9kZW5hIDwtIG5hLm9taXQobW9kZW5hKQoKaW5mZWN0aW9ucy5nYW0gPC0gZ2FtKERhaWx5IH4gcyhUZW1wZXJhdHVyZS4xKSArIHMoSHVtaWRpdHkuMSkgKyBzKHBtMjUuMSkgKyBzKG5waS4xLCBrPTQpLCBmYW1pbHk9cG9pc3NvbihsaW5rPSJsb2ciKSwgZGF0YT1tb2RlbmEpCnBsb3QoaW5mZWN0aW9ucy5nYW0sc2NhbGU9MCxzZT0yLCBzaGFkZT1UUlVFLHBhZ2VzPTEpCnRpdGxlKG1haW4gPSBpKQp9CgpgYGAKCkZsb3JlbmNlOgpgYGB7cn0KZm9yKGkgaW4gMDpudW1faXRlcikgewpudW1fbGFnIDwtIGkKZmxvcmVuY2UgPC0gcmVhZC5jc3YoImZsb3JlbmNlX2Z1bGwuY3N2IikgI1JlYWQgaW4gZmlsZQojUHVzaCB0ZW1wZXJhdHVyZSZodW1pZGl0eSBkYXRhIGJhY2sgY2xpbWF0ZV9sYWcgZGF5cyBieSBhZGRpbmcgZW1wdHkgZGF0YSBpbiBmcm9udAojcmVwKE5BLCBuKSBjcmVhdGVzIGEgdmVjdG9yIG9mIGxlbmd0aCBuIGZpbGxlZCB3aXRoIE5BCmZsb3JlbmNlJFRlbXBlcmF0dXJlLjEgPC0gYyhyZXAoTkEsIGNsaW1hdGVfbGFnKSwgIGZsb3JlbmNlJFRlbXBlcmF0dXJlWygxKToobnJvdyhmbG9yZW5jZSktY2xpbWF0ZV9sYWcpXSkgCmZsb3JlbmNlJEh1bWlkaXR5LjEgPC0gYyhyZXAoTkEsIGNsaW1hdGVfbGFnKSwgZmxvcmVuY2UkSHVtaWRpdHlbKDEpOihucm93KGZsb3JlbmNlKS1jbGltYXRlX2xhZyldKQpmbG9yZW5jZSRwbTI1LjEgPC0gYyhyZXAoTkEsIGNsaW1hdGVfbGFnKSwgZmxvcmVuY2UkcG0yNVsoMSk6KG5yb3coZmxvcmVuY2UpLWNsaW1hdGVfbGFnKV0pCgojcHVzaCBucGkgYmFjayBudW1fbGFnIGRheXMKZmxvcmVuY2UkbnBpLjEgPC0gYyhyZXAoTkEsIG51bV9sYWcpLCBmbG9yZW5jZSRucGlbKDEpOihucm93KGZsb3JlbmNlKS1udW1fbGFnKV0pCiNPbWl0IHRoZSBuZXcgTkFzCmZsb3JlbmNlIDwtIG5hLm9taXQoZmxvcmVuY2UpCgppbmZlY3Rpb25zLmdhbSA8LSBnYW0oRGFpbHkgfiBzKFRlbXBlcmF0dXJlLjEpICsgcyhIdW1pZGl0eS4xKSArIHMocG0yNS4xKSArIHMobnBpLjEsIGs9NCksIGZhbWlseT1wb2lzc29uKGxpbms9ImxvZyIpLCBkYXRhPWZsb3JlbmNlKQpwbG90KGluZmVjdGlvbnMuZ2FtLHNjYWxlPTAsc2U9Miwgc2hhZGU9VFJVRSxwYWdlcz0xKQp0aXRsZShtYWluID0gaSkKfQpgYGAKCgpSb21lOgpgYGB7cn0KZm9yKGkgaW4gMDpudW1faXRlcikgewpudW1fbGFnIDwtIGkKcm9tZSA8LSByZWFkLmNzdigicm9tZV9mdWxsLmNzdiIpICNSZWFkIGluIGZpbGUKI1B1c2ggdGVtcGVyYXR1cmUmaHVtaWRpdHkgZGF0YSBiYWNrIGNsaW1hdGVfbGFnIGRheXMgYnkgYWRkaW5nIGVtcHR5IGRhdGEgaW4gZnJvbnQKI3JlcChOQSwgbikgY3JlYXRlcyBhIHZlY3RvciBvZiBsZW5ndGggbiBmaWxsZWQgd2l0aCBOQQpyb21lJFRlbXBlcmF0dXJlLjEgPC0gYyhyZXAoTkEsIGNsaW1hdGVfbGFnKSwgIHJvbWUkVGVtcGVyYXR1cmVbKDEpOihucm93KHJvbWUpLWNsaW1hdGVfbGFnKV0pIApyb21lJEh1bWlkaXR5LjEgPC0gYyhyZXAoTkEsIGNsaW1hdGVfbGFnKSwgcm9tZSRIdW1pZGl0eVsoMSk6KG5yb3cocm9tZSktY2xpbWF0ZV9sYWcpXSkKcm9tZSRwbTI1LjEgPC0gYyhyZXAoTkEsIGNsaW1hdGVfbGFnKSwgcm9tZSRwbTI1WygxKToobnJvdyhyb21lKS1jbGltYXRlX2xhZyldKQoKI3B1c2ggbnBpIGJhY2sgbnVtX2xhZyBkYXlzCnJvbWUkbnBpLjEgPC0gYyhyZXAoTkEsIG51bV9sYWcpLCByb21lJG5waVsoMSk6KG5yb3cocm9tZSktbnVtX2xhZyldKQojT21pdCB0aGUgbmV3IE5Bcwpyb21lIDwtIG5hLm9taXQocm9tZSkKCmluZmVjdGlvbnMuZ2FtIDwtIGdhbShEYWlseSB+IHMoVGVtcGVyYXR1cmUuMSkgKyBzKEh1bWlkaXR5LjEpICsgcyhwbTI1LjEpK3MobnBpLjEsIGs9NCksIGZhbWlseT1wb2lzc29uKGxpbms9ImxvZyIpLCBkYXRhPXJvbWUpCnBsb3QoaW5mZWN0aW9ucy5nYW0sc2NhbGU9MCxzZT0yLCBzaGFkZT1UUlVFLHBhZ2VzPTEpCnRpdGxlKG1haW4gPSBpKQp9CmBgYAo=